import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt


def myfft_function(signal, begin, end):
    N = end - begin  # 采样点数
    sample_freq = 20  # 采样频率 20 Hz, 大于两倍的最高频率
    sample_interval = 1 / sample_freq  # 采样间隔s
    signal_len = N * sample_interval  # 信号长度
    t = np.arange(0, signal_len, sample_interval)

    fft_data = fft(signal)

    # 计算幅值
    # 在python的计算方式中，fft结果的直接取模和真实信号的幅值不一样。
    # 对于非直流量的频率，直接取模幅值会扩大N/2倍， 所以需要除了N乘以2。
    # 对于直流量的频率(0Hz)，直接取模幅值会扩大N倍，所以需要除了N。
    fft_amp0 = np.array(np.abs(fft_data) / N * 2)  # 用于计算双边谱
    fft_amp0[0] = 0.5 * fft_amp0[0]
    N_2 = int(N / 2)
    fft_amp = fft_amp0[0:N_2]  # 单边谱

    # 计算相位
    phase_spectrum = np.angle(fft_data, deg=True)  # deg=True表示返回的是角度值，否则是弧度值
    #print(phase_spectrum)
    # # 将相位从-180~180转换到0~255
    n180 = np.full(N, 180)  # 创建含有N个数的数组，值都为180
    print(n180)
    phase_spectrum = (phase_spectrum + n180) * 0.7
    #print(phase_spectrum)
    fft_pha = phase_spectrum[0:N_2]

    # 计算频谱的频率轴
    list1 = np.array(range(0, int(N / 2)))
    freq = sample_freq * list1 / N  # 单边谱的频率轴

    return freq, fft_amp, fft_pha, t


# 主函数
N = 1000                       # 采样点数
sample_freq = 20                 # 采样频率
sample_interval = 1/sample_freq   # 采样间隔
signal_len = N*sample_interval    # 信号长度
t1 = np.arange(0, signal_len, sample_interval)

# 生成采样信号
freq_origin = 0.05
origin1 = 1 * np.sin(2*np.pi*freq_origin*t1)  # 采集的信号
origin2 = 1 * np.sin(2*np.pi*freq_origin*t1 + 1)  # 采集的信号

freq1, amp1, pha1, t = myfft_function(origin1, 0, 1000)
freq2, amp2, pha2, t = myfft_function(origin2, 0, 1000)



# 绘制结果
plt.figure()
# 原信号
plt.subplot(221)
plt.plot(t, origin1, 'y')
plt.plot(t, origin2, 'g')
plt.title(' Original signal')
plt.xlabel('t (s)')
plt.ylabel(' Amplitude ')

# 幅度信号
plt.subplot(222)
plt.plot(freq1, amp1, 'y')
plt.plot(freq2, amp2, 'g')
plt.title(' FuPinTu')
# plt.ylim(0, 250)
plt.xlim(0, 1)
plt.xlabel('frequency  (Hz)')
plt.ylabel(' amp ')

# 相位信号
plt.subplot(223)
plt.plot(freq1, pha1, 'y')
plt.plot(freq2, pha2, 'g')
plt.plot(freq_origin, np.interp(freq_origin, freq1, pha1), 'r*')  # 使用np.interp(freq_origin, freq1, pha1)可以获得指定横坐标的纵坐标值
plt.plot(freq_origin, np.interp(freq_origin, freq2, pha2), 'r*')
plt.title(' XiangPinTu')
# plt.ylim(0, 250)
plt.xlim(0, 1)
plt.xlabel('frequency  (Hz)')
plt.ylabel(' pha ')

# 相位差
plt.subplot(224)
plt.plot(freq1, abs(pha2-pha1), 'b')
plt.plot(freq_origin, np.interp(freq_origin, freq1, abs(pha2-pha1)), 'r*')
plt.title(' XiangWeiCha')
# plt.ylim(0, 250)
plt.xlim(0, 1)
plt.xlabel('frequency  (Hz)')
plt.ylabel(' pha_dif ')

plt.show()